Beyond Strain Release: Delocalization-Enabled Organic Reactivity

The release of strain energy is a fundamental driving force for organic reactions. However, absolute strain energy alone is an insufficient predictor of reactivity, evidenced by the similar ring strain but disparate reactivity of cyclopropanes and cyclobutanes. In this work, we demonstrate that electronic delocalization is a key factor that operates alongside strain release to boost, or even dominate, reactivity. This delocalization principle extends across a wide range of molecules containing three-membered rings such as epoxides, aziridines, and propellanes and also applies to strain-driven cycloaddition reactions. Our findings lead to a “rule of thumb” for the accurate prediction of activation barriers in such systems, which can be easily applied to reactions involving many of the strained building blocks commonly encountered in organic synthesis, medicinal chemistry, polymer science, and bioconjugation. Given the significance of electronic delocalization in organic chemistry, for example in aromatic π-systems and hyperconjugation, we anticipate that this concept will serve as a versatile tool to understand and predict organic reactivity.


S1 Computational methods
All energies and thermodynamic quantities reported in this work were obtained using ORCA (v.4.2.1). 1 Minima and transition states (TSs) were initially identified using autodE (v.1.0.0b3), 2 with low energy conformers located using the ETKDGv2 algorithm implemented in RDKit v. 2019.03.4, 3 and optimised using GFN2-xTB implemented in xTB (v 6.2.2) 4 followed by PBE0-D3BJ/def2-TZVP//PBE0-D3BJ/def2-SVP in ORCA (v.4.2.1). 1,5,6Anionic reactions were run using the GBSA 7 / CPCM 8 solvent models for THF in xTB / ORCA, respectively.Geometries and energy were then refined in ORCA at the [DLPNO-CCSD(T)/def2-QZVPP (TightPNO)//B2PLYP-D3BJ/def2-TZVP] level of theory (CH 3 reactions) or [SMD(THF)-DLPNO-CCSD(T)/ma-def2-QZVPP (TightPNO)// SMD(THF)-B2PLYP-D3BJ/def2-TZVP (ma-def2-TZVP on N)] level of theory (NH - 2 reactions). 5,6,9,10ll calculations used the resolution of the identity approximation (RIJCOSX), 11 with the appropriate auxiliary basis sets. 12'Tight' optimisation criteria (10 −8 Ha tolerance for SCF, 10 −6 Ha tolerance for optimisation step) were employed along with Grid6 / GridX6, corresponding to a Lebedev-590 angular grid, and a radial integral accuracy (IntAcc) of 5.34.Stationary points for the model systems were characterised through calculation of the Hessian.Minima were characterised by the absence of imaginary frequencies, and TSs by the presence of a single imaginary mode.Grimme's quasi-RRHO method 13 was used to calculate entropic corrections to obtain free energies at 298.15 K as implemented in the Python package otherm. 14For reactions calculated in the gas phase, a 1 atm standard state was employed.For reactions in implicit solvent, a 1 atm to 1 M standard state correction was applied by adding RT ln(24.5)= 1.89 kcal mol −1 to the calculated free energy of each species.NBO occupation numbers were calculated using the NBO program (v. 7.0), and ELF descriptors were calculated with Multiwfn (v.3.6). 15ll data processing was carried out using the Scikit-learn package in Python 3.7, 16 and MLR plots were generated with Matplotlib. 17A Python script to generate plots is included as part of the Supporting Information.Individual figures can be generated interactively, or to plot all figures using the terminal run: for i in {1..21}; do echo $i | python mlr_models.py-v; done Enthalpies were chosen for a direct comparison with strain energies, which are commonly reported instead of Gibbs free energies.Trends in enthalpy and Gibbs free energy were found to be in excellent agreement for all reactions studied here.Values of 2 − N occ were found to be in good agreement with an alternative density-based delocalization parameter, Dσ D 0 σ , which was evaluated at the bond critical point (Fig. S3). 18

S2 Linear free energy and geometric relationships Activation barrier prediction
The Marcus equation for a chemical reaction is given by eq S1, where ∆E ‡ is the activation barrier for a given reaction and ∆E ‡ int is the intrinsic activation barrier in the absence of a driving force (∆E r = 0).
For two similar reactions with equal driving force, according to eq S1 any difference in their activation barriers will be described by the difference in their intrinsic activation barriers, ∆∆E ‡ int : Combining equations S1 and S2, we arrive at a variant of the Marcus model that accounts for the effect of varying both the reaction driving force (through ∆E r ) and the intrinsic activation barrier (through ∆E ‡ int ) on the observed activation energy (eq S3), relative to a reference intrinsic reaction barrier ∆E ‡ int (0), where ∆E ‡ int = ∆E ‡ int (0) + ∆∆E ‡ int .The sensitivity of ∆E ‡ toward variation in the driving force and intrinsic activation barrier are given by equations S4a and S4b, respectively.
While the value of ∆∆E ‡ int is a priori unknown, we can replace it with a calculated parameter that captures the physical origin of the change in ∆E ‡ int for the range of systems of interest.For example, in the context of small ring reactivity, we propose that a more delocalized bond will be associated with a lower intrinsic activation barrier; in other words, it will be inherently easier to break a more delocalized bond.We can formulate a linear free energy relationship (eq S5) based on this hypothesis in which we use the parameterization ∆∆E ‡ int = κχ, where χ is a variable that captures bond delocalization (vide infra) and κ is a proportionality constant that gives κχ the units of energy.The sensitivity constants defined in equations S4a and S4b are replaced by fitting parameters α and β, respectively, where κ has been absorbed into the β parameter.An increase in driving force is expected to decrease ∆E ‡ (α > 0), and an increase in bond delocalization is also predicted to decrease ∆E ‡ (β > 0).These parameters may be found by MLR using a Bell-Evans-Polanyi-type approximation, through which ∆E r and ∆∆E ‡ int are assumed to be uncorrelated.The quality of this assumption can be assessed through deviation of α and β from 1 2 and κ, respectively.
Transition state geometry prediction We may derive an equation to predict the extension of a breaking bond from its equilibrium value to that found at the TS as follows: The equations of two parabolas, E 1 and E 2 , that describe the diabatic potential energy surfaces of a reaction are given by where k is the force constant associated with the stretch of each bond, r is the length of the breaking bond, r ′ is the hypothetical length at which the breaking bond is fully cleaved, and ∆E r is the reaction driving force.The x coordinate of the intersection point of the two parabolas (∆r ‡ ), representing the position of the TS, is given by Using a Marcus-type approach, we can define 1 2 r ′ as the intrinsic TS bond extension (∆r ‡ int ), for a symmetrical TS in which ∆E r = 0 such that eq S7 becomes This equation represents a Bell-Evans-Polanyi-like equation for the prediction of TS geometric parameters using knowledge of only the reaction driving force, the force constant associated with the breaking bond, and the intrinsic bond extension.
To enable examination of the effect of varying ∆r ‡ int on ∆r ‡ , we may use the substitution ∆r ‡ int = ∆r ‡ int (0) + ∆∆r ‡ int , resulting in eq S9, where ∆r ‡ int (0) is the intrinsic TS bond extension for ∆E r = 0 and ∆∆r ‡ int = 0. Partial derivatives with respect to ∆E r (fixed ∆r ‡ int ) and ∆r ‡ int (fixed ∆E r ) give the sensitivity of ∆r ‡ toward changes in driving force (eq S10a) and intrinsic bond extension (eq S10b), respectively.
Since ∆∆r ‡ int is a priori unknown, we can replace it with a calculated parameter, in the same way as for the prediction of TS barriers in eq S5.For example, to investigate the role of bond delocalization on TS geometry, we can formulate a linear geometric relationship, first using the parameterization ∆∆r ‡ int = λχ, where χ has the same meaning as before, and λ is a proportionality constant that gives λχ the units of distance.Second, the sensitivity constants defined in equations S10a and S10b may be replaced by fitting parameters γ and δ, respectively, where λ has been absorbed into the δ parameter.An increase in driving force is expected to decrease ∆r ‡ (γ > 0), and an increase in bond delocalization is also predicted to decrease ∆r ‡ (δ > 0).As before, we may find these parameters using MLR, where we again make a Bell-Evans-Polanyi-type approximation through which ∆E r and ∆∆r ‡ int are assumed to be uncorrelated.The quality of this assumption can be assessed through deviation of γ and δ from (2k∆r ‡ int (0)) −1 and λ, respectively.

S4 Heterocycle ring-opening reactivity
Based on the experimental data, we can roughly estimate the relative reaction rates for the two reactions at a given temperature (k rel = k D ′ (T )/k E ′ (T )) using the Eyring equation, if we assume that the experimental conditions are identical except for changes in temperature, that the two mechanisms are identical, and that the reactions occur at identical rates at the two different temperatures used in the reactions (i.e., k D ′ (T )/k E ′ (T ′ ) = 1).
Using the latter condition, the relationship between the free energy barriers of the two reactions and the two reaction temperatures is where R is the gas constant.Using T ′ /T = 1.2, we obtain ∆G ‡ E ′ = 1.2∆G ‡ D ′ +0.1 kcal mol −1 .Based on the reported reaction time of 24 h for D ′ at 298 K, we can estimate a value of ∆G ‡ D ′ ≈ 24 kcal mol −1 (t 1/2 = 12 h).The resulting estimate for ∆G ‡ E ′ is then 29 kcal mol −1 -therefore the free energy barrier is ≈ 5 kcal mol −1 lower for D ′ than E ′ , a difference of only 1 kcal mol −1 from the rule of thumb prediction.Despite completely neglecting entropic effects and avoiding any electronic structure calculations or experiments, the rule of thumb gives a good estimate of the expected reactivity difference based only on tabulated data and visual inspection.Table S9: Differences in thermodynamic quantities (kcal mol −1 ) for the cycloaddition between methyl azide and a range of alkynes, at the B2PLYP-D3BJ/def2-TZVP level.

Figure S4 :
Figure S4: MLR plot (CH 3 + hydrocarbon) for the prediction of ∆H ‡ from ∆H r and E HOM O using ∆H ‡ = ∆H ‡ int + α∆H r + βE HOM O , where α and β are optimised coefficients.The blue dashed line denotes perfect correlation.

Figure S5 :
Figure S5: MLR plot (CH 3 + hydrocarbon) for the prediction of ∆H ‡ from ∆H r and E LU M O using ∆H ‡ = ∆H ‡ int + α∆H r + βE LU M O , where α and β are optimised coefficients.The blue dashed line denotes perfect correlation.

Figure S6 :
Figure S6: MLR plot (CH 3 + hydrocarbon) for the prediction of ∆H ‡ from ∆H r and ∆E HOM O−LU M O using ∆H ‡ = ∆H ‡ int + α∆H r + β∆E HOM O−LU M O , where α and β are optimised coefficients.The blue dashed line denotes perfect correlation.

Figure S9 : 2 Figure S10 :
Figure S9: Anionic and radical additions to a set of 18 heterocyclic molecules (Fig. S10, and the error in activation energy (E a ) predicted by Marcus theory vs 2 − N occ .Reaction data from refs.[19, 20].

Figure S11 :
Figure S11: Set of strain release energies (SREs, kcal mol −1 ) and 2 − N occ values (e) per bond type for a range of mono-, bi-and tricyclic ring systems, cyclic alkynes and alkenes.

Table S1 :
Calculated and predicted activation enthalpies (kcal mol −1 ) for the addition of a methyl radical to each of the molecules in the H12 set.Enthalpies calculated at the DLPNO-CCSD(T)/def2-QZVPP (TightPNO) level, with thermal corrections from the B2PLYP-D3BJ/def2-TZVP level.

Table S4 :
Differences in thermodynamic quantities (kcal mol −1 ) for the addition of a methyl radical to each of the molecules in the H12 set, optimised at the B2PLYP-D3BJ/def2-TZVP level, with single point energies calculated at the DLPNO- CCSD(T)/def2-QZVPP (TightPNO) level.∆H and ∆G at the DLPNO-CCSD(T) level calculated using thermal corrections from the B2PLYP level.